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Abstract 

In this paper, we numerically investigate an inverse problem of recovering the potential term in 
a fractional Sturm-Liouville problem from one spectrum. The qualitative behaviors of the eigenval- 
ues and eigenfunctions are discussed, and numerical reconstructions of the potential with a Newton 
method from finite spectral data are presented. Surprisingly, it allows very satisfactory reconstruc- 
tions for both smooth and discontinuous potentials, provided that the order a € (1, 2) of fractional 
derivative is sufficiently away from 2. 

1 Introduction 

We consider the Sturm-Liouville problem (SLP) for the fractional differential equation 

{- Dgu + qu = \u < X < 1, 
u{0) = u{l) = 0, 

where the (left-sided) Caputo fractional derivative Dg of order a G (1, 2) is defined by 

where r(-) refers to the standard Gamma function. We shall normalize the eigenfunction u by u'(0) = 1. 
As the exponent a tends to 2, problem (1) reduces to the classical SLP [If), Thm. 2.1, pp. 92]. The 
fractional SLP (1) is of immense interest for several reasons. 

First, it arises naturally in the analysis of (spatially) fractional wave equations when applying Fourier 
transform. Fractional wave equations are often used to faithfully capture dynamical behaviors of amor- 
phous materials, e.g., polymer and porous media; see the comprehensive reviews [17, 22] and references 
therein. Hence, a better understanding of the SLP (1) would shed valuable insight into the underlying 
physics of these phenological models. 

Second, diffusion is one of the most prominent transport mechanisms found in nature. At a microscopic 
level, it is the result of the random motion of individual particles, and the use of the Laplace operator 
in the canonical model rests on a Gaussian process assumption on the random motion. Over the last 
twenty years a wide body of the literature has shown that anomalous diffusion in which the mean square 
variance grows faster (superdiffusion) or slower (subdiffusion) than in a Gaussian process offers a superior 
fit to experimental data (see the reviews [29, 2(), 2, 14, !■'>])■ This is particularly true in materials with 
memory, e.g., viscoelastic materials. This anomalous diffusion can be in either the spatial or temporal 
variables and is typically reflected in a fractional order derivative. In either case, separation of variables 
leads to a (fractional order) ordinary differential equation with a set of constants, i.e., the eigenvalues. 

Third, it serves as a natural departure from the classical SLP (a ~ 2), for which there is a wealth of 
profound theoretical results for both the forward and inverse problems (see [4, Chap. 3] for an overview). 
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It is clearly interesting to see how these results might be translated to the new scenario both from a 
mathematical perspective and for what this tells us about the underlying physics. 

Suppose that we wish to solve the inverse problem of recovering the potential q(x) from one spectrum 
{Afc}. It is well known that for the classical SLP (a = 2), this is in general insufficient. One spectrum 
completely determines the potential only if additional a priori information is available, for example, if it 
is known to be symmetric about the midpoint of the interval or is given on one half of the interval and 
has to be determined on the other half. In case of a general potential and fixed boundary conditions at 
the endpoints one requires additional information in the form of a second sequence. This can be a second 
spectrum arising from different boundary conditions, specifying both the eigenfunction and its derivative 
at an endpoint, or giving the energy associated with each frequency - the so-called norming constants (see 
[4, Chap. 3]). However there are many situations where this additional knowledge cannot be obtained, 
and one natural question is whether the same requirements also hold for the fractional case, (1). 

A number of efficient reconstruction techniques, including an iterative method involving an associated 
hyperbolic equation [24], variational methods [21], and linearization methods using an iterated Newton 
method based on the coefficient-to-data map have been developed for the classical case. 

The main result of this paper is strong evidence that the fractional SLP (1) has very different properties 
and extensive numerical results lead to the conclusion that one spectrum completely determines a general 
potential q for 1 < a < 2, which contrasts sharply with the classical case of a = 2. In this paper we will 
also deal with reconstruction methods and follow the Newton scheme as the method of choice. 

We note that despite the extensive literature on the forward problem on fractional ordinary/partial 
differential equations (see [8, 19, 10, 11, 3, 25] for an incomplete list) and their diverse physical and 
engineering applications [IG], the research on relevant inverse problems remains very scarce. Recently, 
several works on unique identifiability and numerical methods for inverse problems for time-fractional 
diffusion/ wave equations [5, 12, 31] have appeared and some unusual phenomena over the classical case 
has been observed. 

The paper is structured as follows. In Section 2 wc discuss results about the Sturm-Liouville theory 
for fractional differential operators. In particular, we develop an asymptotic formula for the spectrum 
and indicate the main properties of the eigenfunctions. It turns out that the spectral values are complex 
and so we have both a real and imaginary parts to the eigenfunctions and it is this feature that allows the 
identifiability from a single, albeit complex spectrum. In Section 3 we numerically investigate the inverse 
SLP with a Newton-type method. Numerical results for both smooth and discontinuous potentials clearly 
illustrate the phenomenon of recovering a general potential for (1) from one spectrum. We conclude the 
paper with some future research problems in Section 4. 

2 Sturm-Liouville theory 

In this section, we discuss qualitative behaviors of the spectrum to the SLP (1). We first discuss the case 
of a zero potential, which unlike the situation with a = 2, is quite nontrivial, and then turn to the case 
of a nonzero potential. 

2.1 Differential operator —Dq 

In this part, we collect known results about the spectra of the fractional differential operators —Dq{1 < 
a <2). We shall use the two-parameter Mittag-Lefher function Ea^p[z) defined by 



The two particular versions of Mittag-Leffier function relevant to problem (1) are Ea^2{z) and Ea.a{z). 
The former occurs in the form x£'q.2(— Ax") as the eigenfunctions of (1) with (7 = 0, while the second 
arises in the kernel of the one-sided Green's function for (1) with q~Q. These functions will also appear 
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later in the Jacobian matrix coming from the hnearized problem. Further, the Dirichlet eigenvalues of 
the operator Dq coincide with zeros of £"0,2(2) [18, Thm. 4]. 

If < a < 2 and € min(7r, a7r)), then the function Ea^fj^z) has the following exponential 
asymptotic expansion (see e.g., [cS, Thm. 1.3-4, pp. 5], [10, eq. (1.8.27), p. 43]): as |z| — >■ 00 



\z-^' exp(.i) - r(/3 - ak) + ^ (p^) ' " ""'^^'^^ " 

(2) 

^ 1 1 / 1 \ 

-ETW3T;n7I + ^bv+T ' /^<|arg(z)|<7r. 



V{ii - ak) 2* 



The fractional SLP (1) was studied earlier [7, 18]. In [7], the existence of a solution to such boundary 
value problem was established, and a certain biorthogonal system was constructed from the eigenfunc- 
tions, see also [8, Chaps. 10-12]. In [18], the aforementioned relation between eigenvalues and zeros of 
Mittag-Leffler function was shown. The distribution of the zeros of the function Ea^p{z) is of independent 
interest [8, Chap. 1] [27, 20]. The next result [8, Lemma 1.4-2, pp. 7] represents one of the main results. 
For the sake of completeness, we sketch the proof. 

Theorem 2.1. For all sufficiently large n (in absolute value), the zeros {zn} of the Mittag-Leffler function 
Ea.2{z) are simple, and have the following asymptotics 



z^ = 2n7ri - (a - 1) (^ln27rlrtl 4- -sign(n) ij + In ^ + d„, 
where the remainder dn is O ■ 

Proof. By taking = 1 in the exponential asymptotics (2), we get 

^-^(^) = ^^""^^'-f(^; + ^(i 

for \z\ — > 00. Hence we have 

1 (2 — a) \z ^ 

Next let C = za and w = C + (a — 1) In^. Then the equation can be rewritten as 



r(2 - a) 

The solutions Wn to the above equation for all sufficiently large n satisfy 

a /I 

Wn = 27m i + In — + O — 

r(2-a) Vl'^l 

or equivalently 

a /I 

C„ + (a - 1) lnC„ = 27rni + In — ^ + O ^ 

r(2-a) VfI 

from which we arrive at the desired assertion 

Q! , X ^ , , Ti" . / N .\ ^ / In jnj 



27rn i + In ■ 



r(2 - a) 



{a - 1) (in 2TT\n\ + ^sign(ri) ij + O 



□ 
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We can improve the estimate in Theorem 2.1 by including further terms in the residual tail in the 
asymptotic (2), N > 2, but with considerable increases in the complexity of the result. 

Clearly, if A € C is an eigenvalue of —Dq, then its conjugate A is also an eigenvalue. Hence we shall 
restrict our discussions to n > 0. As a direct consequence of Theorem 2.1, there are only finitely many 
real eigenvalues for any a < 2. The existence of real eigenvalues is only guaranteed for a sufficiently close 
to 2 [20]. Further, asymptotically, the eigenvalues {A„} are distributed as 



-A„)' 



27rn i + In ■ 



a 



{a 



1) I In 2?™ + - 1 



r(2-a) 

Hence, we immediately get the following corollary. 

Corollary 2.1. Asymptotically, the magnitude |A„| and phase arg(A„) of the n^^ eigenvalue A„ are given 
by 



|A„| - ((27rn + (1 - a)f )2 + ((1 - a) ln2nn + In j^^f) " ^ (2™)", 
27rn+(l-a)f (2 - a)7r 



arg(A„) ~ TT — aatan 



{a — 1) In27rr7, — In 



r(2-Q) 



(3) 




(a) magnitude |An 
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(b) phase arg(A„) 



Figure 1: Numerical verification of the asymptotic formula (3): the solid and dashed 
lines represent the true and predicted values, respectively, for a — 1.1, 4/3, 3/2 and 7/4. 



We show in Figure 1 the magnitude and phase of true eigenvalues and their predictions by the 
asymptotic formula (3). The true eigenvalues are calculated by a quasi-Newton method (cf. Appendix 
B), and the values are precise in the first six digits. The magnitudes can be accurately predicted by 
formula (3), except the first few eigenvalues. Hence, it might allow extracting the exponent a directly 
from the sequence of eigenvalues. The approximation of the phase is not as good, cf. Figure 1(b), and 
it is accurate only for very large n, especially for a values close to 2. This is attributed to both crude 
approximation in deriving the asymptotic formula (3) and the slow convergence of the function atan to 
^. This is in stark contrast with the classical case a = 2, where the asymptotic formula is very accurate 
even for relatively small n [4]. The first ten eigenvalues are shown in Figure 2 for a = 1.10 and 1.75. We 
observe that, at least for a not close to unity, the asymptotic formula (3) is a poor approximation for 
small spectral numbers. So one cannot make much use of this from finite spectral data consisting of the 
lowest frequencies and it is precisely such data that is most easily measured in practice. 
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Figure 2: The first ten eigenvalues for a = 1.10 (left) and a = 1.75 (right). Here the 
dashed line represents the asymptotic with angle (1 — ^)7r. 



The eigenfunctions m„ for the operator — -Dq given by w„ = xEa^2{—^nX°')- These functions are 
reminiscent of sinusoids (of. Figure 3), the eigenfunctions for a = 2, but significantly attenuated close to 
X = \. The degree of attenuation strongly depends on fractional order a. It is observed that for a complex 
eigenvalue, the number of interior zeros of the real and imaginary parts of the respective eigenfunction 
Un always differs by one, and that (either real or imaginary part) for two consecutive eigenvalues differs 
by two. However, the number of zeros for the consecutive real eigenvalues differ only by one (see the last 
row of Figure 3) just as for the classical SLP. 

In the classical inverse SLP almost all reconstruction algorithms make use, directly or indirectly, of 
the fact that the eigenvalues are simple and, moreover, for potentials adjusted to have mean zero, there is 
an interval condition: there is exactly one (Dirichlet) eigenvalue in each interval (7i7r, {n + l)7r) [4, Chap. 
3]. In fact a similar result can be shown for the function Ea^p with a = 2 and 1 < /3 < 3 [b, Thm. 1.4-2]. 
However, if 1 < a < 2 then this extremely useful property is either no longer ensured, or indeed does not 
hold. The asymptotic formula shows that all sufficiently large eigenvalues are simple, but it remains an 
open question if this is true in general - even for the case q = Q\ 

In the classical case, the eigenfunctions have the asymptotic form u„ {x) ~ sin( VA^x) / ^/X^ {cos{\/X^x) 
in the case of non-Dirichlet conditions at a; = 0) [ !], and these are very good approximations for moder- 
ately sized smooth potentials. Thus distinguishing eigenfunctions as well as eigenvalues is straightforward. 
Indeed, this clear correspondence also carries over to the potential recovery: if we arc missing information 
about the m*^ eigenvalue, then we are unable to obtain information about the m^^ Fourier mode of the 
potential [9] but the other modes are essentially unaffected. 

In the case a < 2, the picture is quite different. It is known that for a > 5/3 there are at least 
two real zeros (and hence eigenvalues) to the function Ea,2{z) [20, Thm. 1]. Careful numerical ex- 
periments indicate that the first real zeros appear around a = 1.6 (more precisely within the interval 
(1.5991152,1.5991153)) and they occur in pairs (there are 4 real eigenvalues for a = 1.75 as Figure 2 
shows) so there appears to be always an even number of these. However these pairs can be quite close to 
each other and so there is no possible interval condition. For example, with a = 1.5991153, the two small- 
est eigenvalues in magnitude are real with approximate values 14.0024 and 14.0150. Further refinement 
indicates these two eigenvalues are genuinely simple. The eigenfunction of the lower eigenvalue has no 
zeros in (0, 1), while the larger has one single zero and so they are linearly independent. However, the zero 
occurs at X 0.9994 and the supremum norm difference of the two is less than 1.834 x 10"'*. Thus for all 
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practical computational purposes, neglecting one of these pairs will have no effect. With a = f .599025, 
there are no real roots and the two smallest eigenvalues in magnitude are the complex conjugate pairs 
14.0062 ± 0.1955 i. As we noted above, the imaginary part of the corresponding eigenfunction has no 
zeros in (0,1), and the real part does indeed have a zero in this interval. However, this is at a point 
ccq € (0.9997, 0.9998) and the maximum value of the real part on (xg, 1) is less that 5.2501 x 10^*^. Thus 
practically speaking, this zero plays no computational role. These examples also highlight the difficulty 
in obtaining analytic results or estimates on the properties of the eigenvalues or eigenfunctions. 

2.2 Differential operator —D^ + q{x) 

In the classical SLP, the presence of the potential q{x) influences the eigenvalues by 

A„(g) = A„(0)+ / q{t)dt + cn. (4) 
Jo 

where the remainders c„ decay to zero as n — > +oo. The sequence {c„} effectively encodes all the 
information about the potential q{x). This formula represents only a qualitative behavior. In practice, 
the decay rate of c„ can differ markedly, dependent on the smoothness of the potential: the smoother is 
the potential, the faster is the decay. Numerically, we can observe similar behavior for fractional SLPs. 




■ 5 10 15 '0 5 10 15 



n n 

Figure 4: The decay of the remainder c„ for a = 3/2. 



To illustrate the point, we compute the eigenvalues for the following two potentials 

-2x, a:e[0,|], 
9i(.t) =20x3(6-^"-^)' -e-3) and q2{x) = \ ^t+^a-', 

0, elsewhere. 

The profiles of the potentials can be found in Figure 5 shown later in §3.2. The differences between the 
eigenvalues for the differential operator —Dq + q{x) and that of —Dq + q{t)dt are shown in Figure 4 
for a = 3/2. The decay of the remainders c„ is much faster for the smooth potential gi, which holds for 
both real and imaginary parts of eigenvalues. Hence, the first few eigenvalues might allow a reasonable 
reconstruction of the potential in the inversion procedures. In contrast, in case of the discontinuous 
potential 52, the remainder c„ still has noticeably large oscillations for n up to 15. These observations 
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are also valid for other a values, and concur with that for the classical SLP [4. Figure 1, pp. 82]. We 
note that the presence of a nonzero potential may render two neighboring real eigenvalues into a pair 
of complex eigenvalues. For example, in case of a = 1.6, there are two real eigenvalues for q = 0, i.e., 
13.4205 and 14.6454, whereas that for the potentials qi and q2 are given by the complex conjugate pair 
14.6843 ± 1. 71971 and 14.2242 ± 1. 79101, respectively. 

3 An inverse Sturm-Liouville problem 

Now we turn to the inverse problem of recovering the potential q{x) in the fractional SLP (1) from (finite) 
spectral data {Xn}n=i, with the help of a simplified Newton method. Numerically we observe that one 
single spectrum can uniquely determine a general potential. 

3.1 A Newton method 

We shall numerically solve the inverse SLP by a Newton type method based on that of [ ''"] for the classical 
case. First, we denote by u{q, A) the solution to the initial value problem 

-DQU + qu = Xu < X < 1, 

u{0) = 0, it'(O) = 1. ^ ' 

Obviously, the number A G C will be an eigenvalue if the solution u{q, A) satisfies u{q, A)(l) = 0. 

The knowledge of only finite spectral data {Xn}n^i or equivalently the vector A^v ~ (Ai, . . . , Ajv)* S 
C^, as is often the case in practice, precludes the determination of an arbitrary potential q. Instead, we 
seek a potential in a finite-dimensional space spanned by the basis {wk}, and we take its dimension to 
be M < A'^, allowing us the possibility to reduce this value for purposes of regularization, i.e., 

M 

<l^\x) ^^qkWk{x). 
fc=l 

Now the problem is to find q^^ E spandw^.}) such that u{q'^^ , A)(l) vanishes, i.e., u satisfies also the 
desired right-side boundary condition. We define a nonlinear mapping F : i— > by 

F„(Ajv, q'') = u{q'', n^l,...,N. 

For the given set A^r of eigenvalues, we attempt to solve the system of nonlinear equations 

F{AN,q'')^0. 

Thus Newton's method emerges as a natural candidate (see [1.3] for the case a = 2) and we can 
proceed as follows: given an initial guess g° G span({wfe}*f^-^), we repeat the iteration 

q"+l=q"-(F'(AAr,<z"))"l(F(A^,q")) n - 0, 1, . . . . (6) 

where the vector q" = (g", . . . , q"^) G refers to the expansion coefficients of the approximate potential 
g" in the basis {wk}, until a certain stopping criterion is reached. The entries of the Jacobian i^'(Ajv, q) 
can be found by solving a set of fractional differential equations as we show in the next lemma. 

Lemma 3.1. Thenk^^ entry of the Jacobian F'{AN,q), i.e., ^^"g^"''^'* ; is given by v{Xm q,Wk){i) , where 
the function w(A„, q, Wk) satisfies 

- D^v + qv ^ XnV ~ Wku{Xn,q) < a; < 1, 

v{0) = 0, w'(0) =0. ^ ^ 
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Therefore, evaluating either the map F(AN,q) or the Jacobian F'{AN,q) incurs solving initial value 
problems of fractional order (cf. (5) and (7)), which can be numerically accomplished by the predictor- 
corrector method in Appendix A. The main computational effort in applying the Newton method (6) lies 
in computing the Jacobian F'{AN,q). We note that for large A values, the initial value problem is very 
stiff, and a very refined mesh may be needed to resolve the problem to a prescribed accuracy. 

Of particular interest is the special case q = and using the respective eigenvalues Ajv,o- Then the 
solution v{Xn.o,0,Wk) to problem (7) can be expressed in closed form using the Green's function of the 
operator - X [10, eq. (4.1.74), p. 232] 



v{X„,o,0,Wk){x) 



With the approximation F'(Ajv,07 0) in place of the Jacobian F'(Ajv,9"), the iteration (6) leads to 
the canonical frozen Newton's method, which is known to converge (locally) linearly. This reduces the 
computational expense enormously since computing the Jacobian F'{An, <?") is online and very expensive, 
whereas F' (An^jO) can be calculated offline. When solving (6), we stack the real and imaginary parts 
of the matrix F'{An,o, 0) and the vector F{An, g") as follows 



5R(F'(Ajv,o,0)) 
3(F'(Aw,o,0)) 



and 



^iF{AN,q")) 
^{F{AN,q'')) 



In this situation the system has real components only which ensures that the Newton update J~ r" is 
always real. Experimentally, the frozen Newton's method converges fairly fast. As is a general rule for 
applying these type of methods to nonlinear inverse problem, the invertibility of the Jacobian is often 
very difficult to establish theoretically, and this is the case here. The latter obstacle is closely related to 
the lack of a uniqueness result for the inverse SLP problem. Consequently, we are not able to establish the 
convergence of the quasi-Newton scheme and, nor can we derive error estimates for the finite-dimensional 
approximations. However, it is conjectured that F'{An.o, 0) for any a e (1, 2) is invertible in view of our 
experimental experiences. 



3.2 Numerical results and discussions 

Now we are ready to present numerical reconstructions of the potentials qi and q2 from finite spectral 
data A^r. For the spectral data (Ajv) generation and for the inversion step, we solve the initial value 
problem (5) with the predictor-corrector algorithm (see Appendix A) with a mesh size h = le-3 and 
h ~ 1.25e-3, respectively, in order to reduce the most obvious form of "inverse crime", but not completely 
due to the use of identical predictor-corrector algorithm in the forward and inverse steps. The eigenvalues 
are found by a quasi-Newton method in Appendix B. The basis set {wk} is fixed at {sinkir}. All the 
computations were performed with MATLAB R. 7.12.0(2011a) on a dual core desktop computer. 



Table 1: Reconstruction error e for the potentials. 



a 


1.02 


1.1 


6/5 


4/3 


3/2 5/3 


7/4 


9/5 


1.85 


qi {N^8) 


4.217e-3 


4.156e-3 


4.071e-3 


4.013e-3 


4.042e-3 4.690e-3 


9.169e-3 


2.576e-2 


1.218e-l 


q2 (N^IO) 


1.983e-l 


1.983e-l 


1.983e-l 


1.984e-l 


1.988e-l 2.026e-l 


3.331e-l 


1.709e0 


3.469e0 



The numerical results are summarized in Table 1 and Figure 5. In the table, the error e of an approx- 
imation q to the potential q is defined by e = \\q — q\\L'^(o,i)- We have experimented the reconstruction 
algorithm with different numbers of basis and eigenvalues. Although the concrete numbers differ slightly, 
the observations drawn from these results are the same. Hence, we shall present only the results with 
N = M = 8 and N = M = 10 for the potentials qi and q2, respectively. We remark that the convergence 
of the frozen Newton method is steady and fast, effective numerical convergence being generally achieved 
within five iterations, cf. Figure 6 for an illustration. 



9 



Table 2: Condition number of the stacked Jacobian matrix J. 



a 


1.02 


1.1 


6/5 


4/3 


3/2 


5/3 


7/4 


9/5 


1.85 


N 


= M = 


5 


5.313e0 


6.393e0 


8.302e0 


1.236el 


2.095el 


2.716e2 


2.573e3 


4.136e3 


1.193e4 


N 


= M = 


8 


8.267e0 


1.054el 


1.479el 


2.439el 


4.608el 


6.109e2 


8.767e3 


3.397e5 


1.114e5 


N 


= M = 


10 


1.049el 


1.399el 


2.082el 


3.743el 


7.978el 


1.043e3 


1.110e4 


3.298e5 


3.711e5 


N 


= M = 


15 


1.656el 


2.411el 


4.037el 


8.586el 


2.299e2 


3.998e3 


8.524e4 


1.345e6 


1.895e7 



For a in the range (1, f], the error e of the recovered potential qi largely stays constant, and then 
it deteriorates as the exponent a gets closer to 2. Some reconstructed profiles arc shown in the left 
coltimn of Figure 5. The reconstructions are in excellent agreement with the true potential for a G (li §]■ 
Hence, one single spectrum allows accurately reconstructing a general potential. For a value very close 
to 2, e.g., a = 1.9, we are unable to obtain a good approximation by even using 15 eigenvalues. This 
is expected since in the limit case a = 2, it is well known that one single spectrum is insufficient 
for completely determine the potential, and only the symmetric part can be recovered [13]. This is 
corroborated by the dramatic increase of the condition number of the Jacobian F'{An^o, 0), cf. Table 2: 
the conditioning of the stacked Jacobian J worsens steadily as the exponent a approaches 2 as expected 
since it must be singular for a = 2. The ill-conditioning is inherently related to lack of information and 
is particularly relevant for noisy data. There arc several possible remedies. The most straightforward 
one is to use a smaller number of basis functions which directly leads to a better conditioned Jacobian 
matrix. This strategy has proved very effective in our experiments. Another very natural approach 
would be to employ regularization techniques, e.g., Tikhonov regularization and truncated singular value 
decomposition. Perhaps surprisingly, these were not able to produce better results and sometimes resulted 
in inferior ones. 

Similar observations can be drawn from the numerical results for the discontinuous potential 52; see 
the right column of Figure 5. The reconstruction in case of a = | (or in fact, any lesser value of a) is 
very close to the best Fourier approximation of the exact potential, and given our basis representation 
approach, it represents an optimal reconstruction. For larger a, and strongly so as a approaches a = 2, 
the higher-order Fourier coefficients needed for faithfully capturing the potential q2 are less well resolved 
and hence the reconstructions arc less accurate. In comparison with the smooth potential, the onset of 
deterioration also occurs for a smaller a value. 

These empirical observations strongly indicate that the eigenvalues of the fractional SLP contain 
more information than that of the standard second-order SLP, especially for a values sufficiently away 
from 2. This might be attributed to the fact that the eigenvalues of the fractional SLP arc genuinely 
complex, rather than real as for the classical SLP. Similar phenomena have been observed previously in 
related problems in second-order SLPs and first-order systems. In [23] the authors demonstrated the 
unique identifiability of the potential from one complex eigenvalue sequence for the second-order SLP 
with a complex impedance boundary condition. Yamamoto [■>()] showed that two sequences of (complex) 
eigenvalues of a certain first-order matrix operator uniquely determine two unknown coefficients in the 
operator. 

Lastly, wc contrast the experimentally observed uniqueness for the fractional SLP with the classical 
a = 2 in the linearized case. It is tempting to attribute the very different behavior of the fractional case 
1 < a < 2 from the classical a = 2 to the "anomalous diffusion" assumption. But the behavior for a = 2 
might be viewed as the unusual circumstance. To this end, we let qi and 52 be two potentials, and {A„} 
and {fJ-n} be respective Dirichlet spectrum, and {u„} and {«„} cigenfunctions. Then 

/ {Q1 - q2)UnVn dx = (A„ - /i„) / W„W„ dx . 
Jo JQ 

If the potentials have the same spectrum, then 

1 

(gi - q2)unVn dx = Vn. 
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Figure 6: Convergence of the frozen Newton method for a = 3/2. Here res and err refer 
to the residual t|F(g, Ajv)|| and the error e, respectively. 



Now in the case a = 2, the cigenfunctions behave asymptotically like sui{^/X^x) + 0(--^) 

and VAn = nir + and similarly {v„} and /i„ [4, Chap. 3]. Meanwhile, the asymptotic expansion 

implies that the mean values of the potentials qi and q2 are identical. Hence, the identity sin^ nnx = 
i(l — cos2n7rx) implies that to leading order, there holds 

/ {Qi ~ 92) cos(2ri,7rx) dx = Vn. 
Jo 

Therefore, to leading order, the even part of the potential difference qi — (72 is zero; that is the Dirichlct 
spectrum determines the even part of the difference, but equally, gives no information on the odd part. 
A refinement of this idea leads to the original uniqueness proof due to Borg [] ] of both the two spectrum 
problem and the case when the potential is known to be symmetric. 

However, this proof and its consequences rely on the fact that the set of squares of the cigenfunc- 
tions {un} only spans effectively "one half" of the space L^{Q, 1), upon ignoring the small errors in the 
procedure, and this is due entirely to the trigonometric identity relating the squares of cigenfunctions 
of order n to cigenfunctions of order 2n. In the case 1 < a < 2, even for g = 0, we have no such 
identity relating to another single eigenfunction; in fact we would expect an expression of the form 
""raC^) = J^'kLi '^kUkix) whcrc none of the complex-valued sequence {a^} vanishes. Looked at in this 
light, the anomalous behavior is in fact exhibited by the case a = 2. 

4 Concluding remarks 

We have presented a first numerical study of the forward and inverse Sturm-Liouville problems for frac- 
tional differential operators. The numerical results indicate a stunning observation: one single spectrum 
can uniquely determine a general potential, and a few (finite) spectral data allow very good reconstruc- 
tions, provided that the fractional order a £ (1, 2) is sufficiently less than 2. 

These experimental observations lead to a number of interesting conjectures for the fractional Sturm- 
Liouville problem: 

(a) All the eigenvalues of the operator —Dq {a E (1, 2)) are simple. 
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(b) The (Dirichlct) eigenvalues of ~Dq + q{x) (1 < a < 2) obey the relation (4) with {c„} S 

(c) For any a G (1,2), one single (Dirichlet) spectrum completely determines the potential q in the 
fractional Sturm-Liouville problem (1). 

(d) While (c) is the optimal result, a weaker version might be more tractable: to show the injectivity 
of the linearized map. 

(e) There are clear extensions to other boundary conditions and also the case of the half line. 

Apart from these theoretical questions, the rigorous analysis of relevant numerical schemes is also 
of immense interest. These include error estimates for the eigenvalue approximations, convergence and 
convergence rates of the quasi-Newton scheme etc. A better understanding of the analytical aspects of 
the fractional Sturm-Liouville problem is crucial for addressing these numerical issues. 
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A Predictor-corrector algorithm 

Here we describe a predictor-corrector algorithm for fractional initial value problems based on the idea in 
the work [(>]. There are several errors in the formulas in [(>], in, e.g. eq. (15). Hence we present necessary 
details here. Consider the initial value problem (with a G (1, 2)) 

( D'^u^ f{x,u{x)) 0<a;<l, 
\u(0)=uo, u{{))^ull\ 
This initial value problem is equivalent to the Volterra integral equation [10, eq. (3.1.41), pp. 141] 



u(x) = uo + + -— / {:x-tr-'f{t,u{t))dt. (8) 
i (a) Jo 

To solve equation (8), we partition the unit interval [0,1] into a uniform mesh T = {xk = kh,k = 
0,1,..., /C}, X S N, with a mesh size h = The predictor-corrector algorithm in [6] extends the 
Adams-Bashforth-Moulton method. The predictor step applies the left rectangle rule (with the kernel 
(.Tfc+i — x)"""^ being the weight) to the integral in (8), i.e.. 



/ {xk+1 - xY^^g{x)dx « — V bj^k+i9{xj), 

Jo a ,_n 



3=0 

where the integrand g is continuous and the weights &j,fc+i sltc given by 

b,.k+i = (fc + 1 - jT - (fc - jT, J - 0, 1, . . . , fc. 

The corrector step applies the trapezoidal rule (again with the kernel {xk+i — x)"~^ being the weight) to 
the integral in (8), i.e., 

/ {xk+1 - x)''-^g{x)dx « - — — — V aj^k+ig{xj), 
Jo (a + Ij" ^ 
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where the weights cij^k+i are given by 

aj,k+i^{ (fc-j + 2)"+i + (fc-jr+'-2(fc-J + l)"+\ if.? = l,...,fc, 
[ 1, ifj = fc + l. 

The predictor-corrector algorithm apphes repeatedly the above two quadrature rules to get the ap- 
proximations u^{xk+i) and respectively by 
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,(1). 



Then the sought-for approximation U}i{xk-i-i) is obtained by re-evaluate the corrector step with u'j^{xk+i) 
in place of u^{xk+i). To further improve the accuracy of the approximation Uh, we adopt a Richardson 
extrapolation step, which reads Uh = ^(4uh — u^). The solver for the initial value problem always 
incorporates this extrapolation step. 



B Quasi-Newton method for eigenvalue problem 

Next we describe a quasi- Newton method for finding an eigenvalue A to the SLP (1). The starting point 
of the method is the following obvious observation: any eigenvalue A to the SLP (1) is one root of the 
nonlinear map from the potential q to u{q, A)(l), where u{q, A) is the solution to the initial value problem 
(5). Clearly, then the solution u[q, A) will also be the respective eigenfunction. 

We shall find the eigenvalues to the SLP (1) by a quasi- Newton method, the secant method. The 
complete algorithm is listed in Algorithm 1 . The stopping criterion at Step 8 can be based on monitoring 
the quasi-Newton update 5X: if the increment \5\\ falls below a given tolerance (which is set to 1.0 x 10~^^ 
in our computation), then the algorithm is terminated. Accurate initial guesses for the SLP (1) with 
zero potential can be obtained from the zeros of Mittag-Leffler function Ea,2{—X), which can be directly 
extracted by visualizing the function i?Q.2(— A) ['2(S]. That for a nonzero potential can be obtained by 
adding the integral q{t)dt to the eigenvalues for zero potential. In practice, including a small imaginary 
part in the initial guesses is beneficial. Each iteration of the algorithm invokes solving one initial value 
problem, which is dominant in the expense. Our numerical experiences indicate that its convergence 
appears always superlinear, and the convergence basin is relatively large. 



References 

[1] G. Borg. Eine Umkehrung der Sturm-Liouvilleschen Eigenwertaufgabe. Acta Math., 76(l):l-96, 
1946. 

[2] J. -P. Bouchaud and A. Georges. Anomalous diffusion in disordered media: Statistical mechanisms, 
models and physical applications. Phys. Rep., 195(4-5):127-293, 1990. 

[3] H. Brunner, L. Ling, and M. Yamamoto. Numerical simulations of 2D fractional subdiffusion prob- 
lems. J. Comput. Phys., 229(18):6613-6622, 2010. 

[4] K. Chadan, D. Colton, L. Paivarinta, and W. Rundell. An Introduction to Inverse Scattering and 
Inverse Spectral Problems. SIAM, Philadelphia, PA, 1997. 

[5] J. Cheng, J. Nakagawa, M. Yamamoto, and T. Yamazaki. Uniqueness in an inverse problem for a 
one-dimensional fractional diffusion equation. Inverse Problems, 25(11):115002, 16, 2009. 



14 



Algorithm 1 Quasi-Newton method for eigenvalue problem (1). 



1: Given two initial guesses A° and A^; 

2: Calculate M(A°,g)(l) and u{X^,q){l) by the predictor-corrector method; 
3: for fc = 1, . . . , i^r do 

4; Compute derivative d\u{X'' ,q){l) by finite difference 



Compute the quasi-Newton update 5X by 

^(A^g)(l) 
aA^/(A^<z)(l)' 



6: Update the eigenvalue by A''+^ = A*' + 6X; 

7: Compute m(A'^+^, by predictor-corrector method; 

8: Check stopping criterion. 

9; end for 

10: output approximate eigenvalue X^ and eigenfunction u{X^,q). 



[6] K. Diethelm, N. J. Ford, and A. D. Freed. A predictor-corrector approach for the numerical solution 
of fractional differential equations. Nonlinear Dynani., 29(l-4):3-22, 2002. 

[7] M. M. Djrbashian. A boundary value problem for a Sturm-Liouville type differential operator of 
fractional order. Izv. Akad. Nauk Armjan. SSR Ser. Mat, 5(2):71-96, 1970. 

[8] M. M. Djrbashian. Harmonic Analysis and Boundary Value Problems in the Complex Domain. 
Birkhauser Verlag, Basel, 1993. 

[9] H. Hochstadt. The inverse Sturm-Liouville problem. Comm. Pure Appl. Math., 26:715-729, 1973. 

[10] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo. Theory and Applications oj Fractional Differential 
Equations. Elsevier, Amsterdam, 2006. 

[11] Y. Lin and C. Xu. Finite difference/spectral approximations for the time-fractional diffusion equa- 
tion. J. Comput. Phys., 225(2):1533-1552, 2007. 

[12] J. J. Liu and M. Yamamoto. A backward problem for the time-fractional diffusion equation. Appl. 
Anal., 89(11):1769-1788, 2010. 

[13] B. D. Lowe, M. Pilant, and W. Rundell. The recovery of potentials from finite spectral data. SIAM 
J. Math. Anal., 23(2):482-504, 1992. 

[14] M. M. Meerschaert, D. A. Benson, and B. Baumer. Multidimensional advcction and fractional 
dispersion. Phys. Rev. E, 59(5):5026-5028, 1999. 

[15] R. Metzler and J. Klafter. The random walk's guide to anomalous diffusion: a fractional dynamics 
approach. Phys. Rep., 339(l):l-77, 2000. 

[16] R. Metzler and J. Klafter. The restaurant at the end of the random walk: recent developments in the 
description of anomalous transport by fractional dynamics. J. Phys. A: Math. Cen., 37(31):R161- 
R208, 2004. 

[17] R. Metzler and T. F. Nonnenmacher. Space- and time-fractional diffusion and wave equations, 
fractional Fokker-Planck equations, and physical motivations. Chem. Phys., 284(l-2):67-90, 2002. 



15 



[18] A. M. Nahuscv. The Sturm-Liouvillc problem for a second order ordinary differential equation with 
fractional derivatives in the lower terms. Dokl. Akad. Nauk SSSR, 234(2) :308-311, 1977. 

[19] I. Podlubny. Fractional Differential Equations. Academic Press, San Diego, CA, 1999. 

[20] A. Y. Popov. On the number of real eigenvalues of a certain boundary-value problem for a second- 
order equation with fractional derivative. J. Math. Sci., 151(l):2726-2740, 2008. 

[21] N. Rohrl. A least-squares functional for solving inverse Sturm-Liouville problems. Inverse Problems, 
21(6):2009-2017, 2005. 

[22] Y. A. Rossikhin and M. V. Shitikova. Application of fractional calculus for dynamic problems of 
solid mechanics: novel trends and recent results. Appl. Mech. Rev., 63(1):010801, 52 pp., 2010. 

[23] W. Rundell and P. Sacks. Numerical technique for the inverse resonance problem. J. Comput. Appl. 
Math., 170(2):337-347, 2004. 

[24] W. Rundell and P. E. Sacks. Reconstruction techniques for classical inverse Sturm-Liouville problems. 
Math. Camp., 58(197):161-183, 1992. 

[25] K. Sakamoto and M. Yamamoto. Initial value/boundary value problems for fractional diffusion- wave 
equations and applications to some inverse problems. J. Math. Anal. Appl, 382(l):426-447, 2011. 

[26] W. R. Schneider and W. Wyss. Fractional diffusion and wave equations. J. Math. Phys., 30(1):134- 
144, 1989. 

[27] A. M. Sedletskii. Asymptotic formulas for zeros of functions of Mittag-Lcfher type. Anal. Math., 
20(2):117-132, 1994. 

[28] H. Seybold and R. Hilfer. Numerical algorithm for calculating the generalized Mittag-Lefher function. 
SIAM J. Numer. Anal, 47(l):69-88, 2008/09. 

[29] M. F. Shlesinger, B. J. West, and J. Klafter. Levy dynamics of enhanced diffusion: Application to 
turbulence. Phys. Rev. Lett., 58(11):1100-1103, 1987. 

[30] M. Yamamoto. Inverse spectral problem for systems of ordinary differential equations of first order. 
I. J. Fac. Sci. Univ. Tokyo Sect. lA Math., 35(3):519-546, 1988. 

[31] G. H. Zheng and T. Wei. A new regularization method for the time fractional inverse advection- 
dispersion problem. SIAM J. Numer. Anal, 49(5):1972-1990, 2011. 



16 



